MICROCOPY  RESOLUTION  TEST  CHART 

NATIONAL  BURLAU  Of  STANDARDS  I%.<  A 


_ UNCLASSIFIED _ 

SECURITY  CLASSIFICATION  or  THIS  PAGE  (Whm  Data  En(J.«<Q  * 

REPORT  DOCUMENTATION  PAGE  befowDc^ple™2Nform 

TV'fcpdRT  number  .  ^  M  as  I2.  govt  accession  wo.  »  recipient's  catalog  number 


REPORT  DOCUMENTATION  PAGE 

1.  report  numBEr  .  -  ^  j  ft  I2.  govt  accession  wo. 

AFOSRTE-  84-0  148 

«■  TITLE  (tnd  Submit) 

FEASIBILITY  STUDIES  OF  OPTICAL  PROCESSING  OF  IMAG1 
BANDWIDTH  COMPRESSION  SCHEMES 


t.  authors; 

B.R.  HUNT 
R.N.  STRICKIAND 

R.A.  SCHOWENGERDT _ 

>.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Un‘*0.  0-$  Arixona. 

fi )  VS*nt  Arisen*.  2S1  Su  / 

TTT  CONTROLLING  OFFICE  NAME  ANO  ADDRESS 

AFOSR/NE 

Bolling  AFB,  DC  20332 


S.  TYPE  OF  REPORT  A  PERIOD  COVERED 

ANNUAL 

•.  PERFORMING  ORG.  REPORT  NUMBER 
».  CONTRACT  OR  GRANT  NUMBERfaJ 

AF0SR-81-01 70 


to.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  A  WORK  UNIT  NUMBERS 


61102F  2305/BI 


II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  12-  REPORT  DATE 

AFOSR/NE  _ MAY  15,  1983 _ 

Bolling  AFB,  DC  20332  number  of  pages 

_ _ _ _ _ 68 _ 

<4.  MONITORING  AGENCY  NAME  A  AOORESSfll  dllltnnt  /ram  Controlling  Of  tic*)  IS.  SECURITY  CLASS,  (of  this  report) 

UNCLASSIFIED 

15a.  OECLASSIFICATION/DOWNGRADING 
SCHEDULE 

16.  DISTRIBUTION  STATEMENT  (of  thta  Report)  ~~  “  ~~ 

APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


20.  ABSTRACT  (Continue  on  rovot&a  & tdo  If  nocooaary  and  Identify  by  block  number) 

'  This  research  focuses  on  these  three  areas: 

a'V  formulation  of  alternative  architectural  concepts  for  image  bandwidth 
compression,  i.e.,  the  formulation  of  components  and  schematic  diagrams  which 
differ  from  conventional  digital  bandwidth  compression  schemes  by  being 
Implemented  by  various  optical  computation  methods j 

g)  simulation  of  optical  processing  concepts  for  image  bandwidth  compression 
DO  1473  COITION  OF  I  NOV  AS  IS  OBSOLETE  UNCLASSIFIED  (OVER) 


(OVER)' 


FILE  CORY 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (-Whan  Dara  Bntocod) 

84  03  06  110 


UNCLASSIFIED 


SECURITY  CLASSIFJCATIOH  OF  THIS  FAOEfWTun  Dpia  fntarad) 

20.  (Cont.)  , 

^system  performance  sensitivity;  0  • 

cj  saturation  of  optical  processing  for  image  bandwidth  compression  until  the 
overall  state  of  optical  methods  in  image  compression  becomes  equal  to  that  of 
digital  image  compression. 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  this  FAOE(T*Ti.n  Data  Bntarad) 


AFOSR-TR-  84-0143 


Grant  No.  AF0SR-81-01 70 
Annual  Report 

“Feasibility  Studies  of  Optical  Processing 
of  Image  Bandwidth  Compression  Scnemes" 

by 

B .  R  .  Hunt 
R.  N.  Strickland 
R.  A.  Schowengerdt 

May  IS,  1983 


RESEARCH  (AFSC) 

app>V'.  '  •”*  has  fc-en  revis'd  and  Is 

Di  st  r  i ,  ...  i  A , ^  ,  'V#°™  1AV  13  >12. 

MATTHEW  J.  KLRI'EH  ‘ 

L hi  Introduction  Chiaf ,  I.ohnical  Info  ration  Division 

Grant  No.  AF0SR-81-01 70  has  an  objective  which  is  well- 
summarized  by  the  Grant  title:  "Feasibility  studies  of  optical 
processing  for  image  bandwidth  compression  schemes."  It  is  the 
intent  of  research  sponsored  under  this  Grant  to  direct  investi¬ 
gation  into  the  following  issues: 

(a)  formulation  of  alternative  archi techtural  concepts 
for  image  bandwidth  compression,  i.e.,  the  formula¬ 
tion  of  components  ana  schematic  diagrams  wnicn 

di ffer  from  conventional  digital  bandwidth  compres¬ 
sion  schemes  by  being  implemented  by  various 
optical  computation  methods; 

(b)  simulation  of  optical  processing  concepts  for  image 
bandwidth  compression,  so  as  to  gain  insight  into 
typical  performance  parameters  and  elements  of  system 
performance  sensitivity; 

(c)  maturation  of  optical  processing  for  image  band¬ 
width  compression  until  the  overall  state  of  optical 
methods  In  image  comoresslon  becomes  equal  to  that 
of  digital  image  compression. 

It  is  the  last  of  these,  item  (c),  which  represents  the 
continuing  strategic  objective  of  the  efforts  being  carried  on 
under  Grant  No.  AF0SR-81-0170.  It  Is  important  to  remember  that 
the  major  attention  given  to  image  bandwidth  compression  has 
been  for  methods  most  conveniently  Implemented  by  digital  compu- 
tions.  As  flexible  and  multipurpose  are  digital  metnods,  there 
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may  always  be  operational  circumstances,  environments,  or  con¬ 
straints  where  the  availability  of  a  different  technology  Is 
important.  However,  with  the  concentration  upon  digital  compu¬ 
tations,  which  has  characterized  most  research  on  bandwidth 
compression,  alternative  methods  In  optics  have  suffered.  Thus, 
the  purpose  of  research  sponsored  under  this  Grant  is  to  serve 
as  a  source  of  alternatives  for  future  concepts  in  bandwidth 
compression,  so  that  the  environment  for  compression  technology 
need  not  be  dominated  by  one  methodology. 


(II.)  Overview  of  the  Report 

The  research  currently  sponsored  under  Grant  No.  AFOSR-81- 
0T70  consists  of  several  distinct  and  separate  activities.  The 
separate  research  efforts  are  unified  by  a  common  theme:  the 
application  of  optical  processing  for  image  bandwidth  compres¬ 
sion.  Within  this  common  theme,  however,  the  separate  researcn 
projects  are  not  completely  related  to  each  other.  There¬ 
fore,  this  report  is  put  together,  literally,  as  a  number  of 
independent  reports.  The  separate  sections  of  the  report,  which 
follow  this  section,  are  intended  to  be  read  separately  ana  in¬ 
dependently  of  any  other  section.  Each  section  has  its  own 
references  and  its  own  figure  labellings,  for  example. 

The  separate  sections  of  this  report,  and  the  research 
problems  dealt  with  in  each  section,  are  summarized  in  the  fol- 
lowi ng: 

(1)  Data  compression  by  mul ti -spectral  staggered  sampling, 
and  data  reconstruction  by  spatial  and  spectral  inter¬ 
polation  (see  Section  (III.)). 

(2)  Data  compression  by  optical  tomography,  with  data 
recons tructi on  by  optical  convolution  and  back  pro¬ 
jection  (see  Section  (IV.)). 

(3)  Adaptive  data  compression  by  spatial  trans formati ons 
to  create  a  spatially  stationary  image  (see  Section 
(V.)). 

(4)  Improvement  of  the  optical  data  compression  method 
known  as  I0PCM  (see  Section  (VI.)). 


) 


Mul ti spectral  Staggered  Sampling  Oata  Compression 


This  section  reviews  progress  on  the  project  to  achieve  data 
compression  of  mul tl spectral  Imagery  by  sensor  undersampling  and 
subsequent  recons tructl on  of  the  original  scene  utilizing  Inter¬ 
band  redundancy  of  edge  information.  The  technique  assumes  that 
features  possessing  high  spatial  frequencies  are  similar  in  all 
spectral  bands  of  a  scene. 

The  sampling  scheme  considered  in  the  present  research 
consists  of  bands  of  a  mul ti spectral  sensor  which  are  under¬ 
sampled  at  staggered  intervals  in  a  scene.  For  example,  the 
spatial  arrangement  of  pixels  in  an  equally-sampled  four-band 
image,  bands  I-IV,  is: 


As  shown  in  Figure  1,  the  pixels  of  this  mosaic  arrangement  may 
be  directly  transmitted  (PCM)  for  moderate  data  compression,  or 
may  be  further  encoded  (eg.  OP  CM )  for  higher  compression  rates. 
At  the  receiver,  OPCM  decoding  (if  needed)  is  followed  by  re¬ 
construction  of  all  pixels  in  each  band.  To  accomplish  this, 
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the  mosaic  data  Is  subsampled  for  each  band  in  turn;  these  intra¬ 
band  pixels  are  then  spatially  interpolated  to  the  size  of  the 
full  mosaic  scene  to  become  initial  estimates  in  an  iterative 
space-vari ant  reconstruction  using  the  spectral  correlation  be¬ 
tween  the  bands.' 

The  iteration  scheme  for  recon s tructi on  proceeds  as  follows. 
In  each  iteration,  the  pixels  within  a  moving  spatial  window  a  re 
used  to  calculate  statistics  needed  in  predictors  of  the  missing 
bands.  At  each  window  location  t.ne  center  sample  plus  tne  pre¬ 
dicted  values  are  then  written  to  corresponding  spatial  coordi¬ 
nates  in  four  output  images.  The  reconstructs  on  is  then  either 
terminated  or  the  reconstructi ons  are  used  as  approximate  images 
in  another  iteration. 


If  E  is  the  general  band  to  be  estimated,  and  S,  F,  G,.  .  . 
are  the  other  bands  sampled  within  the  window,  then  an  estimate 
of  E  could  be  formed  from  a  linear  prediction  of  all  bands, 

a 

i.e.,  E  ■  r  +  adS  +  bdF  +  cdG  ♦  .  .  wheredS  *  S  -  J,  etc. 

A  MMSE  cri  teri  on,  such  that  <(E-E)  >  is  a  mi  ni  mum,  res  ul  ts  in  the 
solution  for  the  coupling  parameters  in  terms  of. local  band 


variances  and  spectral  covariances: 

r  .2  .  .  c 


<uSdF> 


<dSdF> 


<ASaG> 

<dFdG> 
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<  AtAr  > 
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The  band  Images  will  be  treated  as  ergodic  random  processes  over 
the  extent  of  the  estimation  window,  so  that  the  ensemble  aver¬ 
ages  are  replaced  by  spatial  averages. 

Simulatl ons 

One-,  two-,  and  three-band  predictors  have  been  tested  in 
the  recons tructi on  of  a  256x25o-pi xel  four-band  LANDSAT  MSS 
image.  3ilinear  interpolations  of  each  band  of  the  mosaic  were 
used  as  initial  estimates  in  the  reconstruction.  Figure  2  is  a 
plot  of  normalized  mean-square-error  (NMSE)  versus  estimator  win¬ 
dow  size  for  the  second  MSS  band.  The  results  are  in  agreement 
with  two  expected  tendencies:  larger  variance  in  the  statistical 
estimates  are  associated  with  decreasing  window  size;  and  higher- 
order  predictors  exhibit  greater  sensitivity  to  these  statistical 
errors.  For  comparison,  also  plotted  in  Figure  2  are  the  results 
when  the  original  images  form  the  initial  estimates.  The  re- 
consturction  error  increases  with  window  size,  reflecting  the 
Increasing  nons tati onari ty  within  the  window.  It  is  seen  that 
all  spectral  regressions  result  in  significantly  smaller  error 
than  does  bilinear  Interpolation. 

Without  direct  compression  of  the  mosaic  data,  spectral 
regression  reconstruct!' on  allows  moderate  data  compression  witn- 
out  compression  hardware  at  the  transmitter.  Figure  3  is  a 
plot  of  NMSE  versus  effective  bit  rate  for  1st-  and  2nd-order 
OP CM  (using  laplaclan  quantl zation)  of  the  LANOSAT  image,  com¬ 
pared  with  spectral  regression  and  bilinear  i nterpol ati on .  A: 
an  effective  bit  rate  of  7  bits/pixel  4-  4  bands  *  1.75  bits/ 
pixel,  spectral  regression  is  numerically  comparable  in 
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accuracy  to  2nd-order  OPCM. 

Further  compression  gains  can  be  made  In  combining 
simple  OPCM  of  the  mosaic  array  with  spectral  regression.  As 
shown  In  Figure  4,  spectral  regression  with  lst-order  mosaic 
OPCM  Is  shown  to  be  more  accurate  at  0.5  bits/pixel  than  direct 
2.nd-order  OPCM  at  1.0  bits/pixel.  This  is  a  significant  per¬ 
formance  gain. 

The  relative  visual  qualities  of  spectral  regression  ana 
OPCM  are  demonstrated  in  Figures  5-3.  Figure  5  is  a  deta-'l 
pnotograpn  of  tne  second  MSS  band  original  scene.  The  optimum 
lst-order  spectral  regression  over  a  5  x  5  window  at  1.75  Pits/ 
pixel  Is  shown  In  Figure  6.  By  comparison,  the  Znd-order  OPCM 
reconstruction  at  2  bits/pixel  is  given  in  Figure  7.  The 
degradation  in  the  regression  case  is  primarily  a  random  pnaslng 
in  the  low-contrast  regions  of  the  scene;  the  degradation  in  the 
OPCM  case  Is  chiefly  granularity  noise.  Finally,  mosaic  OPCM 
with  spectral  regression  at  0.5  bits/pixel  is  displayed  in 
Figure  8. 

Current  research  is  examining  the  effects  of  scene  spatial- 
frequency  bandwidth  and  unequal  band-sampling-rates  on  spectral 
regression  performance. 

Interim  results  of  this  research  have  been  presented  at 
the  1982  annual  meeting  of  the  Optical  Society  of  America;  a 
paper  will  also  be  presented  at  the  August  meeting  of  the 
Society  of  Photo-Optical  Instrumentation  Engineers  in  San  Diego, 
Cali  forni a. 
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Figure  6.  One-Band  Spectral  Regression 
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Figure  7.  OPCM,  2  bits/pixel 


Figure  8.  Mosaic  OPCM  with  spectral  regression,  0.5  bits/pixe' 


Optical  Implementation 

The  staggered  sampling  scheme  may  be  Implemented  either 
by  a  single  mosaic  focal  plane  detector  array  or  by  a  set  of 
staggered  linear  arrays  onto  which  the  scene  Is  scanned  by 
mirrors.  Further,  the  elements  may  either  be  of  dissimilar 
construction  and  spectral  sensitivity,  or  may  consist  of  identi¬ 
cal  wide-band  detectors  to  which  are  bonded  color  filter  arrays 
(CFAs).  Solid-state  color  (RG3)  cameras  have  been  constructed 
which  use  a  set  of  adjoining  linear  CCDs,  each  array  oeing 
sensitive  to  only  one  spectral  band.  The  arrays  are  deployed 
in  alternating  color  sensitivity,  with  elements  of  each  stria 
offset  from  the  elements  of  neighboring  strips  [lj.  Alter¬ 
natively,  dye-deposition  CFAs  have  been  bonded  to  sensor  chips 
[2,3,4].  The  CFAs  are  rectilinear  masks  with  repeating  pattern 
G  3  C2];  the  mask  Is  thus  staggered  by  spectral  element  but 
without  sensor  gaps  as  was  the  case  In  [1].  Configurations 
up  to  484x384  elements  with  34  um  x  20  urn  element  size  have  been 
reported  [4].  However,  little  attempt  has  been  made  to  estimate 
Imagery  In  unsampled  bands  from  data  in  a  sampled  band;  in  [3], 
edge  information  in  the  high  resolution  G  band  was  added  ai  recti 


to  the  low  resolution  R  and  3  bands. 


References 


T.  S.  Ochi,  S.  Yamanaka,  Y.  Kanoh,  and  T.  Nishimura.  "A  Device 
Structure  and  Spatial  Spectrum  Spectrum  for  Checker-Pattern 
CCD  Color  Camera",  IEEE  Trans.  Electron  Devices  ,  Vol.  ED-25 
'  pp.  261-256,  1978. 

2.  8.  Bayer,  "Color  Imaging  Array,"  U.  S.  Patent  3971065  (1976) 

3.  P.  Dillon,  0.  Lewis,  and  F.  Kaspar.  "Color  Imaging  System 
Using  a  Single  CCO  Area  Array,"  IEEE  Trans.  Electron  Devices 
Vol .  ED-25,  pp.  102-107,  1973. 

A.  N .  Koike,  I .  Takemoto,  R.  Satoh,  S.  Hanamura,  S.  Nagaha  ra , 
and  M.  Kubo,  "MOS  Area  Sensor:  Part  I  Design  Consi derati on 
and  Performance  of  an  n-p-n  Structure  484  x  384  Element 
Color  MOS  Imager,"  IEEE  Trans.  Electron  Devices,  Vol.  ED-27, 
pp.  1676-1681,  1980. 


I 


(IV.)  Tomography  In  Image  Data  Compression 

Tomography  Is  a  technique  In  which  a  two-dimensional  Image 
may  be  specified  by  a  number  of  one-dimensional  projections  of 
the  Image  made  at  different  angles.  A  projection  is  obtained 
by  Integrating  the  data  in  one  direction  across  the  image, 
analogous  to  sweeping  a  two-dimensional  dust  pattern  on  a  floor 
into  a  line  with  a  wide  broom.  The  original  image,  or  an  ap¬ 
proximation  to  it,  may  be  recons tructed  from  its  projections  in 
a  number  of  ways.  The  most  obvious  way  is  to  "back-project'  tne 
data  in  each  projection  across  the  image  plane  at  the  angle  of 
the  original  projection  (the  broom  analogy  in  reverse  -  sweeping 
the  dust  line  back  over  a  “sticky'*  floor),  producing  two-dimen¬ 
sional  linear  “smears"  which  are  summed  to  form  a  two-dimensional 
image. 

A  study  of  image  data  compression  through  tomography  is 
being  undertaken.  This  project  was  proposed  to  digitally  simu¬ 
late  a  possible  optical  implementation,  in  which  projections 
of  a  real  Image  are  obtained  with  a  rotating,  cylindrical  lens. 
The  additional  observation  Is  that  projections  are  generally 
slowly  varying  from  one  to  the  next,  suggesting  that  such  in¬ 
herent  data  redundancy  should  allow  further  data  comp  ’•es  s  i  on  . 
Basic  Concepts 

It  is  well  known  that  the  one-dimensional  Fourier  transform 
of  a  tomographic  projection  of  an  image  corresponds  to  the  func¬ 
tion  along  a  radial  “slice"  through  the  two-dimensional  Fourier 
transform  of  the  image  [1].  This  allows  another  method  of 
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Image  recons tructl on ,  In  Fourier  space,  and  gives  additional 
Insight  into  the  working  of  tomography.  A  short  proof  of  this 
“correspondence"  theorem  follows. 

Along  the  horizontal  axis,  a  projection  is  defined  by  the 


equati on : 


v*> 


f(x,y )dx 


where  f(x,y)  is  a  two-dimensional  image  function  in  x  and  y,  f„ 
(y)  is  one  of  its  projections  at  angle  9with  9*0  in  this  case. 
The  one-dimensional  Fourier  transform  of  (1)  is: 


f9(v> 


fQ  (y ) . exp ( -2t1  y v ) dy 


where  exp  is  the  exponential  function,  i  is  /-I  and  again  with 
the  angle  subscript  9*0  in  this  case. 

On  the  other  hand,  the  two-dimensional  Fourier  transform  of 


the  image  function  is: 


F  (u 


•*>  ’ll 


f  (x,y)  .exp(-2TTi  (xu  +  yv))  dxdy. 


Comparing  (2)  and  (3),  it  is  obvious  that: 

F0(v)  *  F  (u,v) 

for  9*0,  that  is  u*Q,  proving  the  theorem  in  this  case. 

If  new  axes  are  chosen  x' ,  y*  in  the  image  plane  and  u',  v‘ 
in  the  Fourier  plane,  each  set  rotated  through  an  angles  from 
the  original  set,  the  rotation  matrix  (old  to  new)  is  given  by: 


COS  9 
•sine 


sine 


cos  9 


which  results  in  the  equality  of  the  new  and  old  kernels  of  the 


integral  in  (3),  thus: 


x '  u 1  ♦  y 1  v 1 


xu  ♦  yv. 
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We  have  that 


P0(R)  •  F  (R  cos 9 1  Rslne) 

as  an  alternate  description  of  the  projection  slice  theorem. 

The  recons turctl on  of  an  image  from  its  set  of  projections 
is  relatively  simple,  but  requires  a  careful  treatment  mathe¬ 
matically.  This  method  we  now  derive  in  the  following  paragraphs, 
since  it  is  important  for  the  purposes  of  the  computer  simula¬ 
tions  which  we  will  be  carrying  out. 

Projection  Geometries  in  Data  Compression 

Using  the  correspondence  theorem  we  can  simulate  tomograpny 
digitally  by  applying  radial  filters  to  a  two-dimensional 
(discrete)  Fourier  transform  of  an  image.  A  digital  image,  held 
as  a  rectangular  array  of  data  (in  x,y)  is  Fourier  transformed, 
using  a  two-dimens  i  ona  T ,  fast  Fourier  transform  (FFT)  algorithm. 
The  result  is  a  rectangular  array  (in  u,v).  A  simple,  radial 
filter  is  used  tc  zero  all  data  not  on  the  desired,  radial  lines. 
An  Inverse  two-dimensional  FFT  results  in  an  image  which  has  been 
effectively  reconstructed  from  the  correspond!' ng  tomogtaphic 
projections.  (See  Figure  1  for  filter  properties.) 

By  varying  the  number  of  radial  "slices"  in  the  filter  the 
effect  of  varying  the  number  of  projections  on  image  quality 
can  be  observed.  A  result  of  the  radial  nature  of  the  slices 
Is  that  image  low  frequencies  are  well  represented,  after 
redundandtly ,  due  to  the  finite  "thickness”  (one  pixel)  of  a 
quantized  slice.  In  contrast,  image  high  frequencies  may  have 


large  regions  in  Fourier  space  deleted,  between  the  filter 
slices.  The  redundancy  at  low  frequencies  corresponds  to  the 
redundancy  previously  noted  by  Hunt  that  projections  are  gen¬ 
erally  slowly  varying  from  one  to  the  next.  This  can  be  re¬ 
stated  by  saying-that  the  higher  frequency  components  of  the 
projections  change  more  rapidly  than  the  lower  frequencies  as 
one  changes  the  angle  of  projection.  Thus,  the  second  imoortant 
procedure  is  to  vary  the  frequency  content  of  (at  least  some  of) 
the  projections. 

3ased  on  Pg(R)  *  r(Rcos6,  Rsin9),  we  derive  the  convolu¬ 
tion-back-projection  method.  Consider  the  inverse  Fourier  trans¬ 
form  of  F(U,V)  in  polar  coordinates.  The  polar  coordinates 
we  will  use  are  expressed  by 


and 


x 

y 

u 

V 


rcos$ 
rsi  n$ 

R  cos8 
R  si n8 


,  r  •  /x2  +  y2 


tan<$ 


.  L 


0  <  r<« 

0  TT 


«(^2  +  V^)sgnV,  -«<R<o 
*  tan’  ^  (  V/U)modT  0<,6<t 


(4) 


(5) 


so  the  differential  area  dUdV  *  |R|dRd6  . 

From  these  definitions,  the  polar  coordinate  Fourier  transform 

of  projection  data  can  be  wri tten  as: 

f(x,y)  *  fn(r,p)  »  f  dS  7  P,(R)  j  R  [  expC  jZ^Rcos  (9-$)]dR  (6) 

p  '  .i  8 
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Let  s  ■  r  cos { 9- $ ) ,  and 

98(s)  -  g0(r  cos(e-t))  - J  P0(R) | R|exp[j2ffRs]  -  PQ ( s )  *q ( s ) 


where 


1 


Pe(s)  *  F  { P ^ ( R ) }  is  the  projection  at  angle 

is  the  kernal  function  . 


q(s)  *  F“ 7  f  j  R i > 


then 


f(x,y)  *  fp(r,$)  * J  g 9 ( r  cos ( 5-p) ) de 


Eq .  (7)  is  the  convolution  (or  filtering)  of  projection  p  (s). 
Eq.(8)  is  the  back-projection  of  Eq .  (7)  to  obtain  the  recon¬ 
structed  image  f  (r,$). 

If  we  assume  the  fp(r,<p)  is  band]  inti  ted  by  a  frequency  3, 

then 

•  a 

q(s)  *  F” 7 { |  R  | }  -  J  |R|eJRsdR  ■  J|R|ej2:rRsdR 

«  2B2  sine  28sir  -  82  sine2  Bs -r  .  (9) 

Eq.  (7)  and  Eq .  (8)  constitute  the  convolution-back-projection 
reconstruct  on  method. 

Digital  Simulations 

The  above  equations  are  for  continuous  data.  For  computer 
simulations  digital  samples  must  be  used,  and  new  equations 
must  be  developed. 

Let  the  given  picture  be  digitized  into  a  grid  of  N2  cells. 
We  assume  that  the  Intensity  of  an  entire  cell  is  concentrated 


(7) 
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at  Its  center  as  shown  In  Fig.  1.  Then  the  pixel  A ( I , J )  at 
location  (x,y)  has  a  corresponding  polar  coordinate  (r,$). 

where 

x  -  (1-0.5)  -  (M/2.0)  .  y  *  (J-Q.5)  -  (M/2.0) 

r  «/x2+y^  ,  tan$  *  £  0^o<2ir  (10) 

Let  the  projection  detector  be  P(N0),  which  consists  of  MO  small 
detectors  and  the  length  of  per  small  detector  is  DL.  We  assume 
the  pixel  A(I,J)  contributes  to  two  small  projection  detectors 
at  most.  Refer  to  Fig.  2,  A  FORTRAN  program  for 
obtaining  the  projection  data  P ( NO )  at  the  projection  angle  e 
is  shown  as  follows: 

DIMENSION  A(N,N)  ,  P(N0) 

NTRF  »  ND/2  +  1 
DO  205  I  «  l.ND 
205  P ( I )  -  0.0 

C**-*  9(in  radian)  is  the  projection  angle. 

PL  1  »  A8S ( cos ( 9 ) ) 

PL2  «  A8S ( s 1 n ( 3 ) ) 

PL  «  AHXKPL1  ,PL2)/0L 
HAPL  »  PL/2.0 
00  215  l  -  I  •  N 
00  215  J  •  1,N 


K 

I 


BE*-? 


c«** 


XPRO  ■  ( r*cos ( 5-9 ) / D L 

(r,§)  is  the  polar  coordinate  of  A(I,J). 

N 1  -  XPRO-HAPL  ♦  0.5 

I F  ( XPRO  -  HAPL  AT.  0.0)  N1  ■  XPRO  -  HAPL  _  0.5 
N2  «  XPRO  +  HAPL  +  0.5 

I F  ( XPRO  +  HAPL  .LT.  0.0)  NT  «  XPRO  -  HAPL  -  0.5 
N3  *  N 1  +  NTRF 
N4  *  N 2  +  NTRF 

C***  NTRF  *  NO/2  +  1  to  transfer  -  N D / 2 <_.N 7  <ND/2  into  i<N3<cND. 

I F C I ABS ( N4  -  N3 )  .  LE .  1  ) GOTO  100 
N5  *  N3  +  1 

P(N5)  «  P(N5)  +  A ( I ,  J  ) 

GOTO  215 

TOO  P  ( N  3 )  •  P  ( N  3 )  +  A ( I , J ) * ( NT  +  0.5  -  (XPRO  -HAPL))/PL 

P(N4)  «  P(N4)  ♦  A(I,J)*(XPR0  +HAPL  -  NT  -  0.5J/PL 
C***  A(I,J)  contributes  linearly  to  small  detectors  N3  and  N4. 
215  CONTINUE 

Based  on  the  above,  we  can  carry-out  the  computer  implementa¬ 
tion  of  the  con  vol  uti  on-back-pro  jecti  on  method. 

Step  1:  For  convolution,  using  Eq.  (7), 

9e(s)  *  P3(s  >  *  q ( s )  -  J  p3(s* )q(s-s* ) ds  *  (7) 

Use  discrete  convolution 

g(na,m4)  •  a  J  p(n' a.mdlqCCn-n1 )a] 

n 

NO/2 

g(na,m<i)  «  a  £  f)(n'a,m<i)n£(r,-n' )a]  (11) 

n'«-H0/2 
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Figure  2.  Pixel  A(I,J)  contributes  linearly  to  two  small 


detectors  N3-  aod  N4 . 


where 


d  •  ^  tn<M-l  M  •  the  number  of  projections 

■  Dl  is  the  length  of  per  small  detector. 

-NO/2  <n<N0/2  ,  -N0/2<n ' <ND/2 

p(n'a,mA)  is  the  projection  data  at  view  angles*  ti. 
q(ea)  is  derived  from  Eq .  (5) 

q ( s )  *  2B2  sine  23STV  -  B2sinc2BST 
by  letting  s  *  li,  i:  integer,  B  *j~(Nyquist  rreq.) 


qUa) 


W 


for  i  *  0 


for  l  *  odd 


for  i  a  even. 


Step  2:  For  back-projection,  using  Eq .  (8), 


f(x.y)  *  fp(r,4>)  -  90( r  cos(9-D))d0 


Using  discrete  integral,  Eq.  (8)  can  be  written: 


Aj  (I,j)  »  f(x,y)  *  fp ( r , s) 


g[r  cos (md- s)  ,m-]  (13/ 


where 


x  «  (1-0.5)  -  N/2.0,  y  -  (J-0.5)  -  N / 2 . 1 
y  •  /x2  +  y2  ,  tan$  ■  £  0<$.<2ir 
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V 


} 


£q.  (13)  Involves  Interpolation  for  estimating  g[r  cos (mA-o ) ,mA] 
from  values  of’Eq.  (11)  g(na,  mA). 

If  na<r  cos (mA-b) <( n+1 )a  , 

Me  estimate  g[r  cos (mA-b) ,mA]  by 

g[r  cos  (mA-b)  ,mA]  ■  i: )(  c  0i.lfT1  -  ~  ^  g(na,mA) 

fl 

d 

Step  3:  Using  Eq .  (11),  £q.  (12),  Eq.  (13),  the  computer 
algorithm  operates  as  follows: 

(1)  Al(I.J)  «  fp(  r,o)  -  0.  |<I,JcN 

(2)  Use  Eq.  (11),  (12)  to  calculate  g(na,  mA) 

(3)  A | ( I , J )  ■  A  1  ( I , J )  +  A . g ( r  cos(mA-q) ,mA) 

where  ■  0£m<M-1  , 

repeat  steps  (2)  and  (3)  M  times,  then  the  final  re¬ 
constructed  image  Aj(I,J)  Is  obtained. 

Progress  to  Date 

Some  qualitative  results  have  been  obtained  with  a  256x256 
element  image  of  a  "baboon".  This  Is  quite  a  good  test  image 
for  the  purpose,  as  it  happens  to  have  a  significantly  flat 
spatial  frequency  spectrum  due  to  the  fine  detail  in  the  baocon's 
fur,  and  includes  areas  of  smooth  grey  values  in  the  eyes  ana 
nose.  The  simple,  radial  filter  has  been  applied  to 
this  image,  varying  the  number  of  radial  slices  and  the  frequency 
cut-offs  of  the  slice.  Surprisi ngly ,  the  image  is  still  just  recog¬ 
nizable  with  only  four  filter  slices  or  projection  planes.  For  a 
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good,  qualitative  result  It  seems  to  be  necessary  to  Include  a 
fairly  large,  central  or  low  frequency  region,  with  a  large 
number  of  high  frequency  "rays'',  forming  a  filter  pattern  in 
Fourier  space  like  a  two-dimensional  porcupine. 

It  is  intended  to  continue  the  experiment  in  a  more  quantita 
ti ve  way,  comparing  image  quality  both  visually  and  by  a  measure 
such  as  image  difference  mean  squares.  In  addition,  some  arti¬ 
ficial  test-pattern  images  will  be  generated,  in  order  to  test 
specific  problems  with  the  method. 

Also,  more  complicated  radial  filters  will  be  tested.  For 
example,  if  a  projection  is  obtained  by  integrating  also  in  a 
circular  direction  (eg.  by  "exposing"  the  projection  sensors 
while  the  projection  angle  is  being  changed),  some  of  the  high 
frequency  information  lost  between  the  radial  filter  slices 
will  be  utilized  (as  averaged  data).  In  the  simulation,  the 
value  at  a  point  on  a  radial  slice  in  the  Fourier  domain  would 
be  replaced  by  the  average  of  values  lying  on  a  short  arc  sub¬ 
tended  by  the  (small)  angle  of  integrated  rotation.  If  all 
values  outside  the  slices  are  then  zeroed  as  before,  the  recon¬ 
struction  is  equivalent  to  rotation  of  the  original  image 
through  a  (small)  angle  during  each  projection  "exposure1.  I”, 
on  the  other  hand,  values  outside  the  slices  in  Fourier  space 
are  replaced  by  the  circular  averages,  the  recons tructi on  is 
equivalent  ot  a  rotation  during  summation  while  constructing 
an  Image  by  back-pro jectl on . 

Finally,  Instead  of  averaging  in  a  circular  direction,  a 
maximizing  function  might  be  used,  thus  accentuating  (though 


in  a  s.ightly  incorrect  position  In  Fourier  space)  any  dominent 
frequency  component  which  happens  not  to  fall  on  a  slice. 
Simulation  Results: 

Based  on  the  above  al gori thm.  for  convolution-back  projection 
results  for  simulating  compression  are  shown  in  Fig. 3,  4,  ana  TaD 


M 

Method 

30 

60 

100 

200 

METH  «  1** 

1.5945 

0.5455 

0.4045 

0.3455 

METH  -  2*** 

1 .419% 

0.553% 

0.457% 

0.415% 

Table  1.  NMSE  #  for  different  projections  numbers  M  and  differ¬ 
ent  kernel  functions. 

* :  NMSE  ■  the  normalized  mean-square-e rror .  [Ref.  1] 


NMSE  * 


x^y  C fr( x ,y )  -  f(x,y)]2 
£  f(x.y)2 
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For  method  1,  the  kernel  function  q(ea)  is 


Q(ea) 


1 
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N4akei 
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for  e  ■  0 

for  e  ■  off 
for  e  ■  even 
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For  method  2,  the  kernel  function  q(ea)  is 


q(ea) 


*2 

TTzazC4e^-l ) 


.  1  *  0 ,  +  1 ,  +  2 ,  .  .  . 

LRef.  3] 
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Figure  3.  The  reconstructed  images  of  different  projection  numbers  using  back 
projection  method  1. 


Figure  4.  The  reconstructed  images  of  dl fferent  projection  numbers  using  back- 
projection  method  2. 
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Optical  Implementations 

One  of  the  most  appealing  aspects  of  the  optical  compres¬ 
sion  scheme  is  how  simply  it  can  be  implemented.  A  one-dimen¬ 
sional  ( cy 1 i ndri ca 1 )  lens  performs  an  integration  of  a  scene 
into  a  projection.  The  image  being  projected  can  be  rotated  by 
any  of  several  simple  optical  components,  such  as  a  Dove  orism. 

The  image  projections  can  be  formed  onto  a  line  detector,  such 
as  a  CCD  array.  Thus,  very  simple  and  inexpensive  .hardware  is 
adequate  for  tne  optical  computation  of  tne  image  projections. 

Equally  simple  hardware  can  be  usao  *‘3r  the  reccns  tructi  on 
of  the  image.  This  nas  not  been  noted  before,  but  is  a  point 
which  we  now  discuss*.  First,  we  note  from  the  above  discus¬ 
sions  that  it  is  necessary  first  to  filter  each  projection  be¬ 
fore  back-projection.  This  being  done,  the  back-projection 
Is  obviously  achieved  by  a  cylindrical  lens,  where  in  this  case 
the  lens  is  projecting  a  line  display  of  the  filtered  projec¬ 
tion  into  a  viewing  space.  The  line  display  can  be  made  by  any 
of  several  methods,  e.g.,  a  line  array  of  L.E.O.’s  or  laser 
diodes.  Probably  the  simplest  Is  a  line  of  intensity  written  onto 
the  face  of  a  CRT.  The  rotation  of  this  display  at  different 
angles  is  carried-out  as  before,  e.g.,  an  image  rctator  prism. 

If  the  projections  and  rotations  are  carried  out  at  a  rata  greater 

than  the  visual  persistence  time  of  the  human  retina,  then  a _ 

*  The  optical  reconsturcti on  is  due  to  Donald  Fraser,  who  re¬ 
ceived  support  from  Grant  no.  AFOSR-31 -0 1 70 ,  while  on  leave  from 
CSIRQ,  Australia. 
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real-time  display  will  be  created.  Since  the  original  tomo¬ 
graphic  compression  system  will  be  operaole  at  video  rates,  the 
reconstruction  system  will  be  operable  at  video  rates  as  well. 

A  simple  video-rate  optical  processing  system  for  data  compres¬ 
sion  can  be  configured  with  existing  off-the-shelf  hardware. 
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We  report  here  two  image  processing  applications  of  sta¬ 
tionary  transform  processing  [1],  (also  described  in  the  1981- 
82  Annual  AFOSR  Report).  The  areas  of  concern  are  image  data 
compression,  and  image  restoration.  In  both  areas,  using  the 
stationary  radiometric  and  geometric  transforms  in  conjunction 
with  nonadaptive  processing  algorithms  proves  to  be  superior 
to  the  application  of  nonadaptive  processing  alone.  In 
section  (V.l)  we  show  results  of  data  compression  *nich  a:o;e/e 
0.6  bits/pixel  for  good  quality  imagery.  Also,  a  scheme  for 
hybrid  optical/digital  implementation  is  proposed.  Finally  in 
section  (V.2),  we  briefly  describe  some  new  work  in  which  we  are 
investigating  the  application  of  stationary  processing  to  spa- 
tl ally-variant  image  restoration. 

(V.l)  Oata  Compression 


The  success  of  image  data  compression  methods  is  closely 
linked  to  the  statistical  behavior  of  the  data.  Conventional 
nonadaptive  predictive  coding  employs  a  global  predictor  to 
describe  the  spatial  correlation  of  the  image.  However,  realis¬ 
tic  Imagery  rarely  exhibits  stationary  statistical  properties. 
Consequently,  the  nonadaptive  approaches  to  predictive,  trans¬ 
form,  and  hybrid  coding  are  far  from  optimum.  Superior  per¬ 
formance  Is  characteri sti c  of  the  more  sophisticated  adaptive 
image  coding  techniques  [2].  An  example  is  adaptive  OPCM,  which 


( 

Is  a  form  of  predictive  coding  where  the  coefficients  of  the 
optimum  predictor  are  tuned  to  'the  local  image  statistics. 

This  adaptation  to  statistical  behavior  leads  to  increased  data 
compression  for  given  image  quality,  or  increased  image  quality 
for  the  equivalent  bit  rates  of  nonadaptive  DPCM.  Similar  bene¬ 
fits  of  adaptive  methods  are  found  in  transform  coding  and 
hybrid  coding. 

We  consider  here  the  reverse  approach  to  the  adaptive  cod¬ 
ing  problem.  Starting  with  a  nons tati ona ry  image,  we  propose  to 
use  a  radi ometri c/geometri c  transform  [1]  which  generates  an 
image  with  approximately  wide-sense  stationary  first  ana  secona 
moment  statistics.  Oue  to  this  stationari ty ,  the  transformed 
image  is  matched  to  nonadaptive  coding  processes.  After  coding 
and  transmission  of  the  transformed  image  data  plus  certain 
transform  coefficients,  the  received  image  is  reverse  transformed. 

Stationary  transform  processing  incorporates  two  separate 
stages:  a  radiometric  transform  which  generates  nearly  station¬ 
ary  mean  and  variance;  followed  by  a  geometric  transform,  or 
warp,  to  give  approximately  stationary  autocorrelation  width. 

In  general,  it  will  be  seen  that  coding  using  stationary  trans¬ 
forms  bears  certain  similarities  to  both  predictive  and  trans¬ 
form  coding.  We  show  that  the  radiometric  transform  acts  in  a 
somewhat  similar  fashion  to  the  adaptive  quantizers  found  in 
other  adaptive  coding  methods.  The  geometric  transform  is 
snown  to  be  a  flexible  and  convenient  means  for  implementing 
variable  spatial  resolution  for  reducing  spatial  redundancy 
within  scenes.  We  consider  in  particular  stationary  transforms 
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applied  to  IDPCM  ( 1 nterpol ati ve  DPCM)  data  compression  [3]. 
I0PCM  is  a  coding  technique  which  incorporates  the  essential 
principles  of  conventional  digital  DPCM  coding;  it  has  been 
proposed  for  Implementation  using  hybrid  digital/optical 
hardware.  In  this  paper,  we  suggest  that  stationary  trans¬ 
forms  are  also  suitable  for  hybrid  di  gi  t  a  1  /or  ti  cal  configura¬ 
tions. 

Applying  Radiometric  Transforms  to  Data  Compression 

?CM  coding  generally  yields  images  of  acceptable  quality 
using  a  minimum  of  5  aits,  or  3 2  gray  levels.  :Jse  of  fewer 
bits  results  in  "false  contouring",  especially  in  low  contrast 
scenes,  since  relatively  large  regions  may  be  assigned  to  a 
single  quantum  level.  Adaptive  quantizers  [4],  [5]  have  been 
developed  which  adaptively  expand  low  contrast  regions  into  the 
full  dynamic  range  of  the  quantizer,  thereby  reducing  the  con¬ 
touring  problem.  Our  scheme  for  achieving  lower  bit  rates 
operates  along  similar  lines.  The  data  compression  scheme  is 
outl 1 ned  in  Fig.  1 . 

The  reduction  in  contouring  effects  using  the  radiometric 
transform  permits  us  to  use  lower  bit  rates.  As  Fig.  1  in¬ 
dicates,  we  also  need  to  transmit  the  values  of  u.,  and  .<  for 

Jl 

each  sub-block.  The  total  bits  can  be  calculated  from  the 
equati on 


8. 


N  n  +  16 


(1) 


total 

where  the  original  image  is  of  size  N  by  N,  n  is  th  number  of 
bits  used  in  quantizing  the  transformed  data,  and  M  is  the 
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sub-block  size.  The  factor  of  16  is  due  to  the  two  bytes  used 


to  code  u^  and  K  per  sub-block.  The  number  of  bits/pixel  is 

B/pi  xel  ■  n  +  ^4-  (2) 


From  Eq.  (2),  the  overhead  bit  rate  needed  to  transmit 
the  transform  coefficients  is  16/Mt;  sub-blocks  of  3  by  8 
pixels  require  0.25  bits/pixel,  for  example.  Larger  sub-blocks 
of  16  by  16  pixels  require  only  0.063  bits/pixel;  however,  the 
image  quality  is  slightly  inferior. 

Figures  2(a)-(f)  illustrate  several  stages  of  the  compres¬ 
sion  scheme  in  Fig.  1.  False  contouring  is  exhibited  in  Fig. 
2(b),  which  represents  normal  3  bit  PCM  coding.  To  overcome 
this  distortion,  mean  and  variance  statistics  are  measured  in 
contiguous  sub-blocks  of  8  by  8  pixels  in  the  "Walter"  image 
of  Fig.  2(a);  these  values  are  then  used  in  the  radiometric 
transform  [1]  In  blockwise  fashion  to  generate  Fig.  2(c).  The 
stationary  value  o$  ■  100  Is  sufficient  to  adaptively  expand  the 
dynamic  range  of  each  sub-block  to  fill  the  available  8  bit 
range.  In  Fig.  2(d),  uniform  8  level  quantization  has  been 
performed  at  levels  16,  48,  80,... ,240.  The  inverse  transform 
[1]  yields  Fig.  2(e).  From  Eq .  2  the  final  bit  rate  is  3.25 
bits/pixel.  Figure  2(f)  is  the  end  result  of  similar  processing, 
except  we  have  used  only  4  quantization  levels  after  transforma¬ 
tion.  The  bit  rate  is  2.25  bits/pixel.  Quantization'  distor¬ 
tion  and  a  block  structure  become  apparent  at  this  bit  rate. 

Oata  compression  by  the  radiometric  transform  method  bears  some 
resemblance  to  the  method  of  block  truncation  coding  [6]. 
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8oth  techniques  require  measurement  of  block  statistics,  fol¬ 
lowed  by  a  quantization  scheme  based  on  these  statistics. 
Applying  Geometric  Transforms  to  Data  Compression 

In  [1]  we  discussed  a  geometric  transform  or  warp  which 
generates  approximately  stationary  autocorrel ati on  length  in 
contiguous  sub-blocks  of  an  image.  This  was  achieved  by  locally 
expanding  (by  nonuniform  interpolation)  sub-blocks  of  low  cor¬ 
relation,  for  example,  edges.  Consider  now  an  alternative  way 
of  generating  nea r-s ta ti ona ri ty  correlation:  compress  or  shrink 
regions  of  high  correlation  until  they  possess  the  same  cor¬ 
relation  factor  as,  say,  the  least  correlated  neighborhood. 
Following  the  procedure  outlined  in  [1],  we  select  3s(^s)  ans 
°s(£H)  t0  be  the  max‘,mura  value  of  o  in  NS  and  EW  directions, 
respectively,  over  all  sub-blocks;  and  apply  compression 
vactors,  F,  to  each  sub-block,  determined  by 

o, 


'NS 


NS 


’s(NS) 


*  ( 0  ^  F ^  1  )  i 


(3) 


with  a  similar  equation  for  F^^. 

The  algorithm  in  [1]  Is  used  to  compute  warp  control  points 
from  the  ^  m  s  * F  E  W  values-  A  scheme  for  data  compression  via 
geometric  trans formati on  is  shown  in  Fig.  3.  This  process 
affords  a  method  of  varying  the  spatial  sampling  interval,  or 
resolution,  in  order  to  compensate  for  high  correlation  or 
spatial  redundancy. 

Each  sub-block  of  the  radiometric  transform  data-compressed 
image  (Fig.  2(e))  Is  represented  in  the  geometric  transform  by 


-  0 


the  following  number  of  pixels 


8/sub-block  ■  M2  Fns  Few  (4) 

The  total  pixels  for  an  image  of  size  N  by  N  is 

Btotal  “  w2a^  N  FNSFEW  (5) 

Assuming  that  a  pair  of  bytes  is  sufficient  to  represent  a  pair 
of  control  points,  the  number  of  overhead  bits  is 

3control  points  *  15  ^  +  ^ 

Ignoring  the  factor  of  1  in  Eq.  (6),  the  final  bit  rate  using 
radiometric  and  geometric  transforms  is,  from  Eqs.  (2), 
and  (6), 
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image  data 
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overhead 
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The  bit  rate  calculation  above  assumes  equal-sized  sub-blocks  in 
the  computation  of  the  two  transforms.  In  practice,  we  can 
usually  use  larger,  hence  fewer,  blocks  for  the  geometric  trans¬ 
form  with  little  effect  on  image  quality.  Block  dimension  is 
more  critical  to  image  quality  in  the  radiometric  transform. 
Substituting  new  sub-block  sizes  MR  and  Mg  in  Eq.  (7)  gives 


B/pixal  « 
1e.  if  Mr  •  Mg,  Eq. 
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all  t: 


fnsfzw 
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(8)  reduces  to  Eq.  (7) 


(3) 
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Following  the  procedure  in  [1],  we  measure  the  values  of 
0^5  and  0£W  at  each  sub-block,  and,  from  Eq.  (3),  compute  the 
compression  ratios  F^  and  Fgy.  In  accordance  with  the  threshold 
limitations  on  p^  and  the  ranges  of  F^  and  F^w  are  3:1  and 

2.6:1,  respectively.  Applying  the  geometrical  transform  based  on 
F^s  and  F^y  [1],  we  obtain  the  waroed  version  of  Fig.  2(e)  in 
Fig.  4(a).  Notice  that  regions  of  low  correlation  remain  rela¬ 
tively  unwarped,  compared  with  regions  of  high  correlation  wnicn 
have  been  undersampled  (along  and  axis)  by  a  factor  approacning 
3:1.  The  inverse  warp  produces  Fig.  4\,b). 

The  terms  in  Eq.  (8)  are:  n  *  3,  N  =  256,  «  16,  =  3, 

and  r *  10^.  Tfie  bit  rate  of  Fig.  4(b)  is  therefore 

1.2  bi ts/pi xel . 

Applying  Stationary  Transforms  to  I  DP  CM  Data  Compression 

Irterpoiati ve  Differential  PCM  (IDPCM)  was  originally  con¬ 
ceived  [3]  as  a  vehicle  for  implementing  conventional  DPCM  using 
hybrid  digital/optical  hardware.  The  method  involves  dividing 
the  image  into  low  and  high  frequency  components,  then  taking 
advantage  of  the  statlsitcal  properties  of  each  component  to 
achieve  data  compression.  Specifically,  the  luw  frequency  image 
is  highly  correlated,  consequently  only  every  other  pixel  on 
every  other  line  is  transmitted,  using  5  bits/pixel  (to  avoid 
false  contouring).  The  high  frequency  component  has  a  small 
absolute  range  of  amplitudes,  and  low  correlation,  hence  every 
pixel  Is  transmitted,  except  at  locations  where  the  low  fre¬ 
quency  data  is  transmitted.  However,  it  has  been  shown  [7] 
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that  a  binary,  or  1  bit,  representation  is  sufficient  to  preserve 
subtle  edge  details  when  added  back  to  the  low  frequency  com¬ 
ponent.  The  bit  rate  is  computed  using 

B/pixel  (IOPCM)  -  j  LF  +  |  HF  (9) 

where  LF  and  HF  are  the  bits/pixel  used  in  the  low  and  high 
frequency  image  components,  respecti vely .  Hence,  using  LF  *  5 
and  HF  *  1 ,  good  quality  imagery  can  be  realized  with  a  data  rate 
of  2  bi ts/pi xel . 

IOPCM  has  two  areas  for  potential  improvement  of  bit 
rates : 

(1)  The  5  bit  coding  of  the  iow  frequency  image  is  inef¬ 
ficient. 

(2)  Except  far  the  subsampl Ing  of  the  low  frequency  image, 
no  attempt  is  made  to  compensate  for  spatial  redundancy  (high 
correlatl on) . 

We  address  these  areas  using  the  radiometric  and  geometric  trans¬ 
forms,  respectively.  Figure  5(a)  Illustrates  our  proposed  in¬ 
clusion  of  stationary  transforms  in  IOPCM.  The  geometric 
transform  Is  first  applied  to  reduce  the  spatial  redundancy 
between  correlated  pixels  In  the  original  image.  We  then  apply 
the  radiometric  transform  to-  "adaptively  quantize"  the  low 
frequency  component.  Fig.  5(b)  shows  the  reconstruction  opera¬ 
tion. 

In  computing  the  final  bit  rate,  we  need  to  consider  four 
separate  contributions: 
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(1)  low  frequency  image 
(ii)  high  frequency  image: 


Ciii)  radiometric  transform 
overhead : 
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(iv)  control  points  overhead:  i  g  |  -L  +  j_  i  -  ig  -r  -  s 

*  LM=  _i  "  V 
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Summing  (i)  -  (iv)  and  dividing  by  N  gives  the  final  bit  rate 


B/pixel  (enhanced  IDPCM) 


(10 


Figures  6(a)  and  6(b)  are  the  results  of  data  compression  by 
enhanced  IDPCM  using  sub-blocks  of  8  by  8  and  16  by  16  pixels, 
respecti vely ,  for  the  radiometric  transform  (the  parameter  MR); 
both  use  16  by  16  blocks  for  deriving  the  geometric  transform 
(the  parameter  M^);  Final  bit  rates  are  0.75  and  0.56  bits/pixel, 
respectively.  The  few  artifacts  which  appear  in  Fig.  6(b)  are 
caused  by  the  larger  sub-block  dimension  of  the  radiometric  trans¬ 
form. 

Hybrid  Ooti cal /Pi gi tal  Implementation  of  Enhanced  IDPCM 

To  Implement  the  IDPCM  technique  described  above  using 
radiometric  and  geometric  transforms,  the  following  sequential 
operations  are  required:’ 
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(A)  A  linearly  subsampled  version  of  the  geometrically  and 
radi ometri cally  transformed  image  is  generated,  quantized,  coded, 
and  transmitted.  (Fig.  5(a)).  This  is  the  low  frequency  chan¬ 
nel.  At  the  same  time,  a  subsampled  version  of  the  geometrically 
transformed  image  that  is  not  radi ometri ca  1  ly  transformed  is 
blurred  (i.e.,  interpolated)  to  produce  the  low  frequency  chan¬ 
nel.  This  blurred  image  is  stored. 

(B)  The  blurred  image  from  the  previous  step  is  subtracted 
from  the  geometri cal ly  transformed  original  image,  producing  the 
geometri cal ly  transformed  high  frequency  image.  This  image  is 
quantized,  coded,  and  transmitted. 

(C)  Reconstruction:  At  the  receiving  end,  the  sub-sampled 
image  from  Step  (A)  Is  decoded;  then  the  inverse  radiometric 
transform  is  applied.  (Fig.  5(b)).  This  subsampled  image  is 
then  blurred  (in  Identical  fashion  to  the  blurring  in  Step  (A)) 
to  generate  the  geometrically  transformed  low  frequency  image. 
Note  that  the  geometric  transform  Is  Implicit  in  the  trans¬ 
mitted,  subsampled  Image,  and  that  the  subsampling  is  in  the 
form  of  a  two-dimensional  linear  array.  The  blurring  is  not 
associated  with  the  geometric  transform,  but  simply  fills  in 
the  missing  pixels  between  the  uniformly  sampled  points.  This 
low  frequency  image  Is  added  to  the  decoded  high  frequency  image, 
and  the  Inverse  geometric  transform  applied  to  the  summed  image 
generates  the  final  Image. 

In  [1]  we  gave  details  of  hybrid  optical/digital  hardware 
for  realizing  block  stationary  transforms.  We  do  not  Include 


i 
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hybrid  data  compression  architecture  here,  since  the  principles 
are  essentially  those  of  reference  [1]. 

(V.2)  Image  Restoration 

We  are  currently  investigating  the  following  problems. 

1.  Spati ally- vari ant  deblurring 

2.  Spati al ly-adapti ve  noise  filtering 

Nonadaptive  noise  filtering  Is  hampered  by  the  inherent  non- 
stationarity  of  image  data.  For  example,  the  lowpass  filter 
required  to  achieve  the  desired  reduction  in  noise  power  may 
cause  intolerable  blurring  o'  edge  details.  Cn  the  other  nana, 
a  filter  which  cuts  off  at  a  higher  frequency  may  not  produce 
enough  noise  smoothing.  What  is  needed,  therefore,  is  a  filter 
which  can  adjust  its  cut-off  frequency  according  to  local  image 
detail  l.e.,  cut-off  being  inversely  proportional  to  correlation. 

We  are  approaching  the  problem  by  using  the  geometric  transform 
to  give  Images  of  approximately  stationary  correlation.  Non¬ 
adaptive  filtering  Is  then  performed  on  the  transformed  images, 
followed  by  Inverse  transformatl on.  (The  radiometric  transform 
Is  not  used  In  this  case,  since  stationary  mean  and  variance 
are  of  no  benefit  to  noise-smoothing.)  Stationarity  in  correla¬ 
tion  length  ensures  that  edges  are  preserved  even  for  high 
degrees  of  noise  smoothing.  A  preliminary  result  of  spatially- 
adaptlve  filtering  is  processed  in  Fig.  7.  The  Improvement  gained  by 
adaptive  filtering  Is  marred  by  striations  caused  by  noise  correlations 
inherent  in  the  filtering  of  geometrically  transformed  images.  .*e  are 
currently  Investigating  ways  of  avoiding  this  effect. 
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C a)  Original 


Pig. 2.  Exaopla  of  data  ccmprassion  using  tie  Radiomatric 
Transfora. 
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(d)  Radiometric  Transform  after  8  level  quantitation 
rig.  2  (continued) 


(£)  Repeat  of  2(d) ,(e)  using  4  levels  (2.25  bits/pixel) 
Fig.  2  (continued) 
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ometric  Transform. 
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(a)  G« omecric  Transfora  of  Fig.  2(e). 

,g.  4.  Example  of  Data  Compression  by  Radicmerric/Geoce 
Transforms . 
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(b)  Result  of  data  compression  using  the 

Radicmetric/Geometric  Transforms  (1.2  bits/pixel) 


Fig. 4.  (continued) 


frequency  channels. 


overhead 


Recons 


7.  Spatially  Adaptive  Noise  Filtering, 
top  left  -  original 
top  right  -  original  pljs  noise 
center  -  adaptively-filtered 
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(VI.)  Improved  IDPCM  Data  Compression 

The  first  technique  which  was  successful  In  demonstrating 
an  optical  architecture  for  Image  data  compression  was  the  method 
referred  to  as  I0PCM  ( Interpo lated  Differential  Pulse  Code  Modula¬ 
tion).  The  essence  of  this  method  Is  to  split  the  Image  Into  two 
components:  a  low-frequency  component  which  is  created  by  optical 
subsampling  and  optical  interpolation  (via  incoherent  convolution), 
and  a  high-frequency  component  created  by  subtracting  the  low- 
frequency  component  from  the  original  image.  The  method  can  be 
implemented  by  simple  optical  masks,  mis-focused  apodized  lenses, 
and  video  subtraction  imaging  [1]. 

The  data  compression  performance  of  IDPCM  is  adequate  but 
not  Impressive.  Indeed,  the  overall  Image  quality  at  a  given 
bit  rata  is  often  reminiscent  of  conventional  DPCM  data  compres¬ 
sion.  For  example,  little  if  any  degradation  in  Image  quality 
is  seen  for  Images  compressed  from  8  bits/pixel  to  3  bits/pixel 
by  both  DPCM  and  IOPCM.  Oropplng  the  bit  rate  to  2  bits/pixel 
leads  to  visible  losses  In  image  information  for  both  DPCM  and 
IDPCM.  An  open  question,  therefore,  is  how  to  achieve  better 
performance  from  the  IDPCM  concept.  One  alternative  is  the 
stationarity  trans formati ons  discussed  in  the  previous  section. 

With  these  methods  image  compression  rates  of  less  than  1  bit/ 
pixel  have  been  achieved  for  IOPCM.  The  stationary  transform 
methods  appear  to  be  somewhat  complex  in  system  archi tecture , 
however,  and  a  desire  exists  to  find  improvements  to  IDPCM  which 
are  less  complicated. 
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An  avenue  of  pursuit  for  Improved  IDPCM  is  In  the  Fourier  . 
coding  technique  which  has  been  used  so  extensively  in  digital 
Image  data  compression.  In  this  method  an  Image  Is  decomposed 
Into  blocks  and  the  2-dlmension  FFT  Is  computed  on  each  block. 

The  Fourier  coefficients  produced  in  each  block  are  then  re¬ 
quantized.  The  key  to  data  compression  by  this  method  is  in  the 
requantization.  Typical  Images  show  a  wide  dynamic  range  in 
the  computed  Fourier  spectrum,  with  the  largest  values  occuring 
near  0.  C.  and  the  smallest  values  near  the  Nyquist  frequency. 
Since  the  dynamic  range  is  so  wide,  it  is  not  necessary  to  carry 
the  same  number  of  bits  of  significance  in  quantizing  all  coef¬ 
ficients  of  the  image  FFT.  For  example,  by  a  suitable  pre-as- 
signed  mask,  which  indicates  the  scale-factor  associated  with 
unit  increments  of  quantization  at  each  spatial  frequency  coef¬ 
ficient,  it  is  possible  to  quantize  coefficients  near  D.  C.  at 
8-10  bits  and  coefficients  near  Nyquist  frequency  at  0-2  bits, 
with  little  overall  error.  Using  these  techniques  in  an  adap¬ 
tive  mode.  It  Is  possible  to  create  quantization  schemes  with 
bit  rates  on  the  average  of  less  than  1  bit  per  pixel. 

The  applicability  of  this  technique  to  IDPCM  has  been  tne 
subject  of  the  current  set  of  investigations.  There  is  substan¬ 
tial  difference  between  IDPCM  coding  of  imagery  and  conventional 
Fourier  coding  of  Imagery,  and  this  difference  leads  to  ques¬ 
tions  as  to  the  applicability  of  Fourier  coding  to  IDPCM.  Soec- 
iflcally,  Fourier  coding  Integrates  the  low  spatial  frequency 
and  high  spatial  frequency  components  Into  one  set  of  quantization 


assignments  for  the  Fourier  frequency  components.  However,  an 
IDPCM  Image  is  split  Into  two  parts,  a  separate  low-  and  high- 
spatial  frequency  set  of  Image  samples. 

Analysis  of  the  IDPCM  technique  Indicates  the  efficiency 
with  which  the  Vow  spatial  frequency  component  is  represented. 

For  example,  consider  an  image  with  8  bits  per  sample  initial 
data  rate.  If  the  subsample  rule  required  for  IDPCM  [1]  is 
chosen  to  be  every  other  pixel  and  line,  then  the  effective  rate 
for  representation  of  low  frequencies  is  8/4  *  2  bits/pixel. 
Further,  it  has  been  shown  that  the  low  frequency  subsamples  need 
not  be  retained  at  the  8-bits/pixel  of  the  original,  and  4  bits 
per  subsample  leads  to  high-quality  results  at  an  average  rate 
of  1  blt/plxel.  Thus,  the  low-spatial -frequency  component  can 
now  be  represented  in  IDPCM  with  rates  of  1  bit/pixel  or  less  [1]. 
This  Indicates  the  area  for  concentrating  upon  Fourier  coding  is 
the  high-frequency  component  of  IDPCM. 

Analysis  of  Fourier  coding  schemes  indicates  that  the  high¬ 
est  spatial  frequencies  In  an  Image  are  representable  with  a 
very  small  number  of  bits.  An  example  can  be  seen  in  Pratt's 
book  [2],  where  a  typical  bit  assignment  mast  for  a  Fourier 
coding  algorithm  averages  less  than  1-bit  per  pixel  for  the 
high-frequencies.  Indeed,  It  Is  aften  possible  to  neglect  the 
highest  spatial  frequencies  In  both  x  and  y  directions.  That 
Is,  a  typical  bit  assignment  mask  may  assign  one-bit  per  Fourier 
coefficient  to  the  x  or  y  Nyqulst  frequency,  but  zero-bits  to 
the  joint  Nyqulst  frequency  of  x  and  y.  (See  Fig.  (VI. 1)). 
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In  Figure  VI. 1  we  see  the  dominance  of  low  frequencies  for  the 
assignment  of  quantization  bits.  In  an  IOPCM  system,  the  low 
frequencies  are  subtracted  out  from  one  Image,  so  that  the  high 
frequency  components  would  represent  only  the  values  that 
would  require  quantization  coding.  A  mask  for  the  Fourier  coding 
of  the  high-frequency  spatial  component  might  appear,  for  example, 
as  in  Figure  VI. 2: 

0  0  0  0  2  1  1  1 

0  0  0  0  2  1  1  1 

0  0  0  0  2  1  1  1 

0  0  0  0  2  1  1  0 

2  2  2  2  2  1  0  0 

111110  0  0 
1  1  1  1  0  0  0  0 

1  1  1  0  0  0  0  0 

The  overall  bit-rate  for  the  high-spatial  frequencies  coded  with 
this  bit  assignment  would  be  much  less  than  one  bit  per  pixel 
(l.e.,  0.677).  Combined  with  an  efficient  sub-sampling  and  re¬ 
quantization  of  the  low  spatial  frequency  components,  a  bit  rate 
of  substanti al ly  less  than  1-blt  per  pixel  should  be  achievable, 
with  acceptable  image  quality. 

The  design  of  the  bit  mask  for  a  Fourier  coding  scheme  can 
not  be  done  by  analytic  model.  The  current  method  for  bit  mask 
by  trial -and-error  studies  over  a  variety  of  images  which  are 
selected  to  be  typical  of  the  class  of  images  to  be  processed 
through  the  data  compression  system.  Thus,  additional  research 
on  this  problem  becomes  one  of  simulation,  experiment,  and  visual 


verification.  Initial  results  Indicate  that: 

(1)  The  Fourier  coding  of  IOPCM  Images  should  produce  images 
of  visual  quality  equal  to  that  of  conventional  Fourier  coding 
of  Images; 

(2)  The  bit-rates  associated  with  a  Fourier  coding  of  IOPCM 
will  be  competitive  with  conventional  Fourier  image  coding. 

Since  IOPCM  was  an  optical  method  for  image  data  compres¬ 
sion,  the  introduction  of  Fourier  coding  leads  to  a  hybrid 
optical /di gi tal  method  for  image  data  compression.  The  reference 
to  a  ‘’hybrid"  method  stems  from  the  fact  that  the  most  straignt- 
forward  means  to  implement  the  Fourier  coding  would  be  by  a 
digital  Fourier  transform  of  the  hi gh-f requency  component  of 
IOPCM. 
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